Results

Transferring Data from Slicer to R

# Load necessary libraries #
library(devtools)
library(SlicerMorphR)
library(rgl)
library(geomorph)
library(readxl)
library(ks)

##############################################
#   Load SlicerMorph Data and Preprocessing  #
##############################################

# Choose the Log File from Slicer
# (Ensure a line break is added at the end of the .log file before loading)
SM.log.file <- file.choose()
SM.log <- parser(SM.log.file, forceLPS = TRUE)

# Extract Landmark Data
Data <- SM.log$LM

# Create geomorph dataframe
datagdf <- geomorph.data.frame(landmarks = Data)  

# Perform Procrustes Superimposition - PCoords stand for Procrustes coordinates
Pcoords <- gpagen(datagdf$landmarks)

# Save Procrustes Coordinates
save(Pcoords, file = "Pcoords.bin")

##############################################
#    Semi-Landmark Sliding (SSL vs HL)       #
##############################################

# Extract Semi-Landmarks Data
Data2 <- SM.log$semiLMs

# SSL: Sliding using **bending energy minimization**
PcoordsBent <- gpagen(A = Data, surfaces = as.numeric(Data2), ProcD = FALSE, print.progress = TRUE)
save(PcoordsBent, file = "PcoordsBent.bin")

# HL: Sliding using **Procrustes distance minimization**
PcoordsProcDist <- gpagen(A = Data, surfaces = as.numeric(Data2), ProcD = TRUE, print.progress = TRUE)
save(PcoordsProcDist, file = "PcoordsProcDist.bin")

# Extract and Save Centroid Size (Csize)
write.table(PcoordsBending$Csize, file = "PcoordsList_SSL.tsv", sep = "\t")
write.table(PcoordsProcDist$Csize, file = "PcoordsList_HL.tsv", sep = "\t")

##############################################
#       Merge with Master Dataset            #
##############################################

# Load Pcoords Lists
PcoordsList_SSL <- read.table(file = 'PcoordsList_SSL.tsv', sep = '\t', header = TRUE)
PcoordsList_HL <- read.table(file = 'PcoordsList_HL.tsv', sep = '\t', header = TRUE)

# Load Master Dataset
Masterfile <- read_excel("path/to/YourMasterfile.xlsx")

# Merge Pcoords with Master Dataset
Merge_SSL <- merge(PcoordsList_SSL, Masterfile, by = "LabID", sort = FALSE)
Merge_HL <- merge(PcoordsList_HL, Masterfile, by = "LabID", sort = FALSE)

# Save Merged Tables
write.csv(Table_SSL, file = "TursiopsSpec_SSL.csv", row.names = FALSE)
write.csv(Table_HL, file = "TursiopsSpec_HL.csv", row.names = FALSE)


##############################################
#   Create and Save Final Analysis Dataset   #
##############################################

# SSL Dataset
spec_SSL <- read.csv("TursiopsSpec_SSL.csv", header = TRUE)
tursiopsBent.dt <- list(gpa.sh = PcoordsBent$coords, cs = PcoordsBent$Csize, spec = spec_SSL)
save(tursiopsBent.dt, file = "tursiopsBent.dt.bin")

# HL Dataset
spec_HL <- read.csv("TursiopsSpec_HL.csv", header = TRUE)
tursiopsProcDist.dt <- list(gpa.sh = PcoordsProcDist$coords, cs = PcoordsProcDist$Csize, spec = spec_HL)
save(tursiopsProcDist.dt, file = "tursiopsProcDist.dt.bin")

##############################################
#         Symmetry Analysis (SSL vs HL)      #
##############################################

## Define Landmark Pairs for Symmetry Analysis

# SSL: Surface Semi-Landmarks - Replace by your own landmarks
lm.pairs_SSL <- matrix(c(
  1:104, 106:109, 111:190, 192:195, 197, 198, 200:217, 219:236, 
  238:287, 290:299, 301, 302, 304:403, 405:434, 436:445, 448:473, 
  475, 476, 479:582, 584:621, 624, 625, 627:640, 642:661, 663:760
), ncol = 2, byrow = TRUE)

# HL: Homologous Landmarks - Replace by your own landmarks
lm.pairs_HL <- matrix(c(1:2, 9:28, 31:56, 59:70, 72:73), ncol = 2, byrow = TRUE)

## Perform Symmetry Analysis

# SSL
asym_SSL <- bilat.symmetry(PcoordsBent$coords, ind = dimnames(PcoordsBent$coords)[[3]], 
                           object.sym = TRUE, land.pairs = lm.pairs_SSL, iter = 9)
symm.sh_SSL <- asym_SSL$symm.shape

# HL
asym_HL <- bilat.symmetry(PcoordsProcDist$coords, ind = dimnames(PcoordsProcDist$coords)[[3]], 
                          object.sym = TRUE, land.pairs = lm.pairs_HL, iter = 9)
symm.sh_HL <- asym_HL$symm.shape

# Save Symmetry Data - sym stands for symmetry
YOURNAME_sym_SSL.dt <- list(gpa.sh = PcoordsBent$coords, cs = PcoordsBent$Csize, 
                                   symm.sh = symm.sh_SSL, spec = spec_SSL)
save(YOURNAME_sym_SSL.dt, file = "YOURNAME_sym_SSL.dt.bin")

YOURNAME_sym_HL.dt <- list(gpa.sh = PcoordsProcDist$coords, cs = PcoordsProcDist$Csize, 
                                  symm.sh = symm.sh_HL, spec = spec_HL)
save(YOURNAME_sym_HL.dt, file = "YOURNAME_sym_HL.dt.bin")

PCA SSL

# Load necessary libraries
library(geomorph)
library(readr)
library(ks)
library(rgl)

# Set working directory (modify as needed)
setwd("INSERT_PATH_HERE")

# Load dataset
load("YourName_sym.dt.bin")


######################################################################
#                          PCA Analysis                              #
######################################################################

# Perform Principal Component Analysis (PCA)
YourNamePCA <- gm.prcomp(YourName_sym.dt$symm.sh)

# Access eigenvalues
eigenvalues <- YourNamePCA$sdev^2

# Calculate proportion of variance explained by each principal component
variance_explained <- eigenvalues / sum(eigenvalues)

# Print proportion of variance explained by each principal component
print(variance_explained)

######################################################################
#                          3D PCA Plot                               #
######################################################################

# Define colors based on Units
group_colors <- rainbow(length(unique(YourName_sym.dt$spec$YourGroup)))

# 3D scatter plot of first three PCs
plot3d(
  YourNamePCA$x[, 1:3], 
  col = group_colors[as.factor(YourName_sym.dt$spec$YourGroup)], 
  size = 20
)

# Add legend
legend3d(
  "topright", 
  legend = unique(as.factor(YourName_sym.dt$spec$YourGroup)), 
  col = group_colors, 
  pch = 10
)

######################################################################
#                        3D Cloud Plot                               #
######################################################################

# Extract PCA scores
Score <- YourNamePCA$x[, 1:3]
Labels <- factor(x = YourNamePCA$YourGroup)

# Kernel density estimation for group classification
Hscv1.Bd <- Hkda(
  x = YourNamePCA$x[, 1:3], 
  x.group = as.factor(YourName_sym.dt$spec$YourGroup), 
  bw = "scv", 
  pre = "sphere"
)

# Perform kernel discriminant analysis
Tt.kda3.Bd <- kda(
  x = YourNamePCA$x[, 1:3], 
  x.group = as.factor(YourName_sym.dt$spec$YourGroup), 
  Hs = Hscv1.Bd
)

# Plot 3D Cloud
plot(
  Tt.kda3.Bd, 
  colors = c("YourColou1","YourColou2","YourColou3","etc.."), 
  drawpoints = TRUE, 
  col.pt = c("YourColou1","YourColou2","YourColou3","etc.."), 
  box = FALSE, 
  display = "rgl"
)

######################################################################
#                       PCA Lollipop Graphs                          #
######################################################################

PC1 <- plotRefToTarget(
  YourNamePCA$shapes$shapes.comp1$min, 
  YourNamePCA$shapes$shapes.comp1$max, 
  method = "vector", 
  mag = 1, 
  gridPars = gridPar(pt.size = 0.25, pt.bg = "red"))

PC2 <- plotRefToTarget(
  YourNamePCA$shapes$shapes.comp2$min, 
  YourNamePCA$shapes$shapes.comp2$max, 
  method = "vector", 
  mag = 1, 
  gridPars = gridPar(pt.size = 0.25, pt.bg = "blue"))

PC3 <- plotRefToTarget(
  YourNamePCA$shapes$shapes.comp3$min, 
  YourNamePCA$shapes$shapes.comp3$max, 
  method = "vector", 
  mag = 1, 
  gridPars = gridPar(pt.size = 0.25, pt.bg = "green"))


######################################################################
#                       PERMANOVA Analysis (PAST)                    #
######################################################################

# Export PCA Scores for further analysis
write.table(YourNamePCA$x, file = "FileName1.tsv", sep = "\t")

A local Image A local Image

A local Image
A local Image